% NEURON_TREE   Export tree as NEURON file.
% (trees package)
%
% [name path] = neuron_tree (intree, name, res, options)
% ------------------------------------------------------
%
% Saves a complete tree in the section based neuron '.hoc' format.
% Alternatively tree can also be stored as '.nrn' file in which each
% segment from the tree graph becomes an independent section in neuron.
%
% Inputs
% ------
% - intree   ::integer: index of tree in trees or structured tree
% - name     ::string:  name of file including the extension ".hoc" or ".nrn"
%     spaces and other weird symbols not allowed!
%     {DEFAULT : open gui fileselect}
% - res      ::number or vector: number of segments per compartment. If
%     vector, one value per compartment, otherwise single value same for
%     all compartments
%     {DEFAULT : ceiling of length in um of compartment}
% - options  ::string:
%     '-s'   : write procedures to collect
%     '-w'   : waitbar
%     '-e'   : include passive electrotonic parameters
%     '->'   : send directly to windows (necessitates -s option)
%     {DEFAULT : ''}
%
% See also load_tree swc_tree start_trees (neu_tree.hoc)
% Uses root_tree cyl_tree dissect_tree ver_tree D
%
% Output
% ------
% - name     ::string:  name of output file;
%     [] no file was selected -> no output
% - path     ::sting: path of the file
%     complete string is therefore: [path name]
%
% Example
% -------
% neuron_tree (sample_tree);
%
% the TREES toolbox: edit, generate, visualise and analyse neuronal trees
% Copyright (C) 2009 - 2017  Hermann Cuntz

function [tname, path] = neuron_tree (intree, tname, res, options)

% trees : contains the tree structures in the trees package
global       trees

if (nargin < 1) || isempty (intree)
    % {DEFAULT tree: last tree in trees cell array}
    intree   = length (trees);
end

ver_tree     (intree); % verify that input is a tree structure

% use full tree for this function
if ~isstruct (intree)
    tree     = trees{intree};
else
    tree     = intree;
end

if (nargin < 4) || isempty (options)
    % {DEFAULT: no option}
    options  = '';
end

% defining a name for the neuron-tree
if (nargin < 2) || isempty (tname)
    [tname, path] = uiputfile ({ ...
        '.hoc', 'export to hoc';...
        '.nrn', 'export to nrn'}, 'Save as', 'tree.hoc');
    if tname  == 0
        tname = [];
        return
    end
else
    path     = '';
end
format       = tname  (end - 3 : end); % input format from extension:
% extract a sensible name from the filename string:
nstart       = unique ([0 (strfind (tname, '/')) (strfind (tname, '\'))]);
name         = tname  (nstart (end) + 1 : end - 4);
name1        = tname  (nstart (end) + 1 : end);
if nstart (end) > 0
    path     = [path tname(1 : nstart (end))];
    tname (1 : nstart (end)) = '';
end
name2        = [path 'run_' name '.hoc']; % show file, with '-s' option
% add a starting node in root to avoid all starting branch point:
tree         = root_tree (tree);

if (nargin < 3) || isempty (res)
    % {DEFAULT: resolution depends on lenght of segment}
    res      = ceil (len_tree (tree));
end

ipar             = ipar_tree (tree); % parent index structure
idpar            = ipar (:, 2);      % index to direct parents
D                = tree.D;           % local diameter values
N                = size (D, 1);      % number of nodes in tree
if isfield       (tree, 'R')
    R            = tree.R;           % region values on nodes in the tree
else
    R            = ones (N, 1);      % add a homogeneous regions field
end
% starting and end points of all branches:
sect             = dissect_tree (tree);
Rsect            = R (sect (:, 2));  % region attribute to sections
uR               = unique (R);       % sorted regions
luR              = length (uR);      % number of regions

if isfield       (tree, 'rnames')
    rnames       = tree.rnames (uR);
    for counter  = 1 : length (uR)
        rnames{counter} = [name '_' rnames{counter}];
    end
else
    if luR       == 1
        rnames   = {name};
    else
        rnames   = cell (1, luR);
        for counter = 1 : luR
            rnames{counter} = [name '_' (num2str (uR (counter)))];
        end
    end
end

nextline         = [(char (13)), (char (10))];
switch           format
    case             '.nrn' % save file in the NEURON .nrn format
        H1           = histc (R, uR); % histogram of regions
        H1 (R (1))   = H1 (R (1)) - 1;
        [~, i2]      = sort (R); % i2 : index of sorted regions...
        iR           = ones (N, 1);
        for counterR = 1 : luR
            if uR (counterR) == R (1)
                iR (1 + find (R (2 : N) == uR (counterR))') =  ...
                    (1 : H1 (counterR))';
            else
                iR (find (R == uR (counterR))') = ...
                    (1 : H1 (counterR))';
            end
        end
        % file-pointer to the neuron-file
        neuron       = fopen ([path tname], 'w');
        % HEADER of the file
        % declaring the regions
        if luR       > 1
            for counterR = 1 : luR
                fwrite (neuron, [ ...
                    'create ' name '_' (num2str (uR (counterR))) ...
                    '[' (num2str (H1 (counterR))) ']', ...
                    nextline], 'char');
            end
        else
            fwrite (neuron, [ ...
                'create ' name '[' (num2str (H1 (counterR))) ']', ...
                nextline], 'char');
        end
        % Connections:
        HW           = waitbar (0, 'writing connections  ...');
        set          (HW, 'Name', '..PLEASE..WAIT..YEAH..');
        fwrite       (neuron, ['proc topolneuron() {', nextline], 'char');
        o            = 0;
        for counter  = 3 : N
            o        = o + 1;
            waitbar  (counter / N, HW);
            if luR   > 1
                fwrite (neuron, [ ...
                    name '_' (num2str (tree.R (idpar (counter)))) ...
                    '['], 'char');
            else
                fwrite (neuron, [ ...
                    name ...
                    '['], 'char');
            end
            fwrite (neuron, [ ...
                (num2str (iR (idpar (counter)) - 1)) ...
                '] connect '], 'char');
            if luR   > 1
                fwrite (neuron, [ ...
                    name '_' (num2str (tree.R (counter))) '['], 'char');
            else
                fwrite (neuron, [ ...
                    name '['], 'char');
            end
            fwrite (neuron, [ ...
                (num2str (iR (counter) - 1)) '] (0), 1',  ...
                nextline], 'char');
            if o == 200
                o = 0;
                fwrite (neuron, [ ...
                    'topolneuron'      (num2str (counter)) '()', ...
                    nextline], 'char');
                fwrite (neuron, [ ...
                    '}', ...
                    nextline], 'char');
                fwrite (neuron, [ ...
                    'proc topolneuron' (num2str (counter)) '() {', ...
                    nextline], 'char');
            end
        end
        close        (HW)
        fwrite       (neuron, ['}', nextline], 'char');
        % Cylinder-Geometry
        fwrite       (neuron, [ ...
            'proc geometry() { local counter', ...
            nextline], 'char');
        for counterR = 1 : luR
            fwrite   (neuron, [ ...
                '   for counter = 0,' (num2str (H1 (counterR) - 1)) ...
                ' {',    nextline], 'char');
            if luR   > 1
                fwrite (neuron, [ ...
                    '      ' name '_' num2str(uR(counterR)) ...
                    '[counter]{', nextline], 'char');
            else
                fwrite (neuron, [ ...
                    '      ' name ...
                    '[counter]{', nextline], 'char');
            end
            fwrite   (neuron, [ ...
                '         pt3dclear()',        nextline], 'char');
            fwrite   (neuron, [ ...
                '         nseg = fscan()',     nextline], 'char');
            fwrite   (neuron, [ ...
                '         pt3dadd(fscan(),fscan(),fscan(),fscan())', ...
                nextline], 'char');
            fwrite   (neuron, [ ...
                '         pt3dadd(fscan(),fscan(),fscan(),fscan())', ...
                nextline], 'char');
            if strfind (options, '-e')
                % passive properties :
                if isfield (tree, 'Ri')
                    fwrite (neuron, [ ...
                        '         insert pas', nextline], 'char');
                    if isfield (tree, 'Ri')
                        fwrite (neuron, [ ...
                            '         Ra = ',    (num2str (tree.Ri)), ...
                            nextline], 'char');
                    end
                    if isfield (tree, 'Gm')
                        fwrite (neuron, [ ...
                            '         g_pas = ', (num2str (tree.Gm)), ...
                            nextline], 'char');
                    end
                    if isfield (tree, 'Cm')
                        fwrite (neuron, [ ...
                            '         cm = ',    (num2str (tree.Cm)), ...
                            nextline], 'char');
                    else
                        fwrite (neuron, [ ...
                            '         cm = 1', ...
                            nextline], 'char');
                    end
                end
                fwrite (neuron, [ ...
                    '         e_pas = 0', nextline], 'char');
            end
            fwrite (neuron, ['      }', nextline], 'char');
            fwrite (neuron, ['   }',    nextline], 'char');
        end
        
        if strfind   (options, '-e')
            % global passive properties if needed:
            fwrite   (neuron, ['forall insert pas', ...
                nextline], 'char');
            fwrite   (neuron, ['forall Ra = ', ...
                (num2str (tree.ri)),         nextline], 'char');
            fwrite   (neuron, ['forall g_pas = ', ...
                (num2str (1 ./ tree.rm)),    nextline], 'char');
            fwrite   (neuron, ['forall cm = ', ...
                (num2str (tree.cm * 1e6)),   nextline], 'char');
            fwrite   (neuron, ['forall e_pas = 0', ...
                nextline], 'char');
        end
        % link to the connections
        fwrite       (neuron, ['topolneuron()', nextline], 'char');
        % footer
        fwrite       (neuron, ['}',             nextline], 'char');
        fwrite       (neuron, ['geometry()',    nextline], 'char');
        % DATA:
        bindex       = 1;
        for counterR = 1 : luR
            if uR (counterR) == R (1)
                bindex = bindex + 1;
            end
            HW       = waitbar (0, [ ...
                'writing cylinders of region' ...
                (num2str (uR (counterR))) ' ...']);
            set      (HW, ...
                'Name',        '..PLEASE..WAIT..YEAH..');
            for counter  = 1 : H1 (counterR)
                waitbar  (counter / H1 (counterR), HW);
                fwrite   (neuron, [ ...
                    (num2str (res (i2 (bindex)))), ' ', ...
                    (num2str (tree.X (idpar (i2 (bindex))))), ' ', ...
                    (num2str (tree.Y (idpar (i2 (bindex))))), ' ', ...
                    (num2str (tree.Z (idpar (i2 (bindex))))), ' ', ...
                    (num2str (tree.D (i2 (bindex)))), ' ', ...
                    (num2str (tree.X (i2 (bindex)))), ' ', ...
                    (num2str (tree.Y (i2 (bindex)))), ' ', ...
                    (num2str (tree.Z (i2 (bindex)))), ' ', ...
                    (num2str (tree.D (i2 (bindex)))), nextline]);
                bindex   = bindex + 1;
            end
            close    (HW);
        end
        fclose       (neuron);
    case             '.hoc' % save file in the NEURON .hoc format
        % distribution of section regions in H1:
        H1           = histc (Rsect, uR);
        % file-pointer to the neuron-file
        neuron       = fopen ([path tname], 'w');
        % HEADER of the file
        fwrite       (neuron, [ ...
            '/*',                nextline], 'char');
        fwrite       (neuron, [ ...
            'This is a CellBuilder-like file written' ...
            ' for the simulator NEURON',        nextline], 'char');
        fwrite       (neuron, [ ...
            'by an automatic procedure "neuron_tree"' ...
            ' part of the TREES package',       nextline], 'char');
        fwrite       (neuron, [ ...
            'in MATLAB',         nextline], 'char');
        fwrite       (neuron, [ ...
            'copyright 2009-2017' ...
            ' Hermann Cuntz',    nextline], 'char');
        fwrite       (neuron, [ ...
            '*/',                nextline], 'char');
        fwrite       (neuron, [ ...
            '',                  nextline], 'char');
        fwrite       (neuron, [ ...
            'proc celldef() {',  nextline], 'char');
        fwrite       (neuron, [ ...
            '  topol()',         nextline], 'char');
        fwrite       (neuron, [ ...
            '  subsets()',       nextline], 'char');
        fwrite       (neuron, [ ...
            '  geom()',          nextline], 'char');
        fwrite       (neuron, [ ...
            '  biophys()',       nextline], 'char');
        fwrite       (neuron, [ ...
            '  geom_nseg()',     nextline], 'char');
        fwrite       (neuron, [ ...
            '}',                 nextline], 'char');
        fwrite       (neuron, [ ...
            '',                  nextline], 'char');
        % declaring the regions
        for counterR = 1 : luR
            fwrite   (neuron, [ ...
                'create ' rnames{counterR} '[' ...
                (num2str (H1 (counterR))) ']', ...
                nextline], 'char');
        end
        fwrite       (neuron, [ ...
            '', nextline], 'char');
        % topology procedure
        fwrite       (neuron, [ ...
            'proc topol_1() {', nextline], 'char');
        countero     = 1;
        counteri     = 1;
        for counter  = 1 : size (sect, 1)
            s        = sect (counter, 1); % start compartment of section
            e        = sect (counter, 2); % end compartment of section
            ipsect   = find (s == sect (:, 2)); % parent section
            ip       = sect (ipsect, 2); % parent index of section
            if ~isempty(ip)
                ie   = find (counter  == find (Rsect == R (e)));
                ipe  = find (ipsect   == find (Rsect == R (ip)));
                fwrite (neuron, [ ...
                    '  connect ', ...
                    rnames{find(uR == R (e))}  ...
                    '[' num2str(ie  - 1) '](0),' ...
                    rnames{find(uR == R (ip))} ...
                    '[' num2str(ipe - 1) '](1)', ...
                    nextline], 'char');
                countero = countero + 1;
                if countero  == 250
                    countero = 1;
                    counteri = counteri + 1;
                    fwrite (neuron, [ ...
                        '}', nextline], 'char');
                    fwrite (neuron, [ ...
                        'proc topol_' ...
                        (num2str (counteri)) '() {', ...
                        nextline], 'char');
                end
            end
        end
        fwrite       (neuron, [ ...
            '}',                               nextline], 'char');
        fwrite       (neuron, [ ...
            'proc topol() {',                  nextline], 'char');
        for counter  = 1 : counteri
            fwrite   (neuron, [ ...
                '  topol_' (num2str (counter)) ...
                '()',                          nextline], 'char');
        end
        fwrite       (neuron, [ ...
            '  basic_shape()',                 nextline], 'char');
        fwrite       (neuron, [ ...
            '}',                               nextline], 'char');
        fwrite       (neuron, [ ...
            '',                                nextline], 'char');
        fwrite       (neuron, [ ...
            'proc shape3d_1() {',              nextline], 'char');
        countero     = 1;
        counteri     = 1;
        for counter  = 1 : size (sect, 1)
            s        = sect (counter, 1); % start compartment of section
            e        = sect (counter, 2); % end compartment of section
            ie       = find (counter == find (Rsect == R (e)));
            fwrite   (neuron, [ ...
                '  ' rnames{find(uR == R(e))} ...
                '[' num2str(ie - 1) '] {pt3dclear()', nextline], 'char');
            indy     = fliplr (ipar (e, 1 : find (ipar (e, :) == s)));
            for counter1 = 1 : length (indy)
                fwrite   (neuron, [ ...
                    '    pt3dadd(', ...
                    (num2str (tree.X (indy (counter1)))),', ', ...
                    (num2str (tree.Y (indy (counter1)))),', ', ...
                    (num2str (tree.Z (indy (counter1)))),', ', ...
                    (num2str (tree.D (indy (counter1)))),')',  ...
                    nextline], 'char');
                countero = countero + 1;
                if countero  == 250
                    countero = 1;
                    counteri = counteri + 1;
                    fwrite (neuron, [ ...
                        '  }',           nextline], 'char');
                    fwrite (neuron, [ ...
                        '}',             nextline], 'char');
                    fwrite (neuron, [ ...
                        'proc shape3d_' num2str(counteri) '() {', ...
                        nextline], 'char');
                    fwrite (neuron, [ ...
                        '  ' (rnames{find(uR == R(e))}) ...
                        '[' (num2str (ie - 1)) '] {', ...
                        nextline], 'char');
                end
            end
            fwrite   (neuron, [ ...
                '  }',              nextline], 'char');
        end
        fwrite       (neuron, [ ...
            '}',                    nextline], 'char');
        fwrite       (neuron, [ ...
            'proc basic_shape() {', nextline], 'char');
        for counter  = 1 : counteri
            fwrite   (neuron, [ ...
                '  shape3d_' (num2str (counter)) '()', ...
                nextline], 'char');
        end
        fwrite       (neuron, [ ...
            '}',                    nextline], 'char');
        fwrite       (neuron, [ ...
            '',                     nextline], 'char');
        fwrite       (neuron, [ ...
            'objref reg_' name, '_all', ...
            nextline], 'char');
        for counter  = 1 : luR
            fwrite   (neuron, [ ...
                'objref reg_' (rnames{counter}), ...
                nextline], 'char');
        end
        fwrite       (neuron, [ ...
            'proc subsets() { local counter', ...
            nextline], 'char');
        fwrite       (neuron, [ ...
            '  reg_' name '_all = new SectionList()', ...
            nextline], 'char');
        for counter  = 1 : luR
            fwrite  (neuron, [ ...
                '  reg_' rnames{counter} ' = new SectionList()', ...
                nextline], 'char');
            fwrite  (neuron, [ ...
                '  for counter = 0, ' ...
                (num2str (H1 (counter) - 1)) ' ' ...
                rnames{counter} '[counter] {', ...
                nextline], 'char');
            fwrite  (neuron, [ ...
                '    reg_' rnames{counter} '.append()', ...
                nextline], 'char');
            fwrite  (neuron, [ ...
                '    reg_', name, '_all.append()', ...
                nextline], 'char');
            fwrite  (neuron, [ ...
                '  }', ...
                nextline], 'char');
        end
        fwrite       (neuron, [ ...
            '}',                   nextline], 'char');
        fwrite       (neuron, [ ...
            'proc geom() {',       nextline], 'char');
        fwrite       (neuron, [ ...
            '}',                   nextline], 'char');
        fwrite       (neuron, [ ...
            'proc geom_nseg() {',  nextline], 'char');
        fwrite       (neuron, [ ...
            '}',                   nextline], 'char');
        fwrite       (neuron, [ ...
            'proc biophys() {',    nextline], 'char');
        fwrite       (neuron, [ ...
            '}',                   nextline], 'char');
        fwrite       (neuron, [ ...
            'access ' rnames{1} ,  nextline], 'char');
        fwrite       (neuron, [ ...
            'celldef()',           nextline], 'char');
        fwrite       (neuron, [ ...
            '',                    nextline], 'char');
        fwrite       (neuron, [ ...
            '',                    nextline], 'char');
        if strfind   (options, '-e')
            fwrite   (neuron, [ ...
                'forsec reg_' name '_all insert pas', ...
                nextline], 'char');
            fwrite   (neuron, [ ...
                'forsec reg_' name '_all g_pas = ' ...
                (num2str (tree.Gm)), ...
                nextline], 'char');
            fwrite   (neuron, [ ...
                'forsec reg_' name '_all Ra = ' ...
                (num2str (tree.Ri)), ...
                nextline], 'char');
            fwrite   (neuron, [ ...
                'forsec reg_' name '_all cm = ' ...
                (num2str (tree.Cm)), ...
                nextline], 'char');
            fwrite   (neuron, [ ...
                'forsec reg_' name '_all e_pas = 0', ...
                nextline], 'char');
            fwrite   (neuron, [ ...
                '',  nextline], 'char');
        end
        fclose       (neuron);
    otherwise
        warning      ('TREES:IO', 'format unknown');
        return
end

if strfind       (options, '-s')
    % file-pointer to the run-file
    neuron       = fopen (name2, 'w');
    fwrite       (neuron, [ ...
        'load_file ("nrngui.hoc")', nextline], 'char');
    fwrite       (neuron, [ ...
        'xopen ("' name1 '")',      nextline], 'char');
    fclose       (neuron);
    if strfind   (options, '->')
        if ispc  % this even calls the file directly (only windows)
            winopen  (name2);
        end
    end
end




